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Abstract—The Macro Basis Functions (MBFs) approach is a 
form of domain-decomposition method applied to radiation and 
scattering problems solved by using integral-equation techniques. 
It enables a systematic reduction of the number of degrees of 
freedom, from that imposed by the discretization of the surfaces 
to that associated with the physical limits of field distributions. 
This paper reviews different variants of this approach, including 
the techniques for determining the MBFs and for fast calculation 
of their interactions. The link with Krylov-subspace iterative 
methods is described, the relationship between the surface of 
subdomains and the number of physical degrees of freedom is 
discussed and multi-level schemes are revisited. Finally, avenues 
for further research are outlined in the Conclusions section of 
this paper. 


Index Terms—macro basis functions, integral equations, 
method of moments, characteristic basis functions, synthetic 
functions 


I. INTRODUCTION 


Efficient and accurate solution of electromagnetic-field inte- 
gral equations has been an important research topic for many 
years. Despite the availability of computers with fast CPUs 
and abundant as well as affordable memory resources, ever- 
increasing demand for solving larger problems still outpaces 
the rapid advances in numerical techniques. The challenges 
faced almost a decade ago were described in a review paper 
[1] and the domain-decomposition approach was introduced 
around the same time frame to solve large problems by using 
the “divide and conquer” approach. One such methodology 
is based on expressing the solutions in the subdomains in 
terms of high-level basis functions that are linear combina- 
tions of a number of pre-computed solutions for the isolated 
subdomains, or for those domains surrounded by relatively 


small neighborhoods. Such a divide and conquer concept 
was already present in earlier works such as [2]-[4] and 
has been developed more systematically by Suter and Mosig 
[5], who introduced the expression “Macro Basis Function” 
(MBF). Quite a few other methods based on aggregation 
of low-level basis functions, such as [6], [7] appeared in 
the computational electromagnetics (CEM) literature almost 
contemporaneously, or soon thereafter. The main attribute of 
the domain decomposition approach is that it enables us to 
handle considerably larger problems, in terms of number of 
Degrees of Freedom (DoFs) than is possible by using the 
conventional Method of Moments (MoM). 

Our objectives in this paper are to review some of the 
earlier works, present the latest developments in this area and 
provide new perspectives on this class of methods. In order to 
facilitate the understanding of the following sections, Section 
II briefly describes what may be viewed as an elementary MBF 
approach, while Section II provides a summary of associated 
methods. Section IV explains how MBFs can be generated, 
while Section V describes different techniques for fast calcula- 
tion of interactions between the MBFs. Following this, Section 
VI reviews the link between MBF approaches and modern 
iterative techniques and Section VII addresses the important 
challenges encountered when attempting to solve multi-scale 
problems. Finally, Section VIII briefly summarizes this work 
and presents some perspectives on the future directions. 


II. ELEMENTARY MBF APPROACH 


This section summarizes what may be regarded as the 
simplest possible MBF approach. For reasons that will be 
apparent later, it may not be the most effective numerical 
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approach, but it will be used to introduce the terminologies and 
notations, while laying the foundations of what is to follow. 
We will assume that the reader is already familiar with the 
Method of Moments (MoM). 

Let us write the original MoM system of equations as Z x = 
v in which Z is the MoM impedance matrix, x is the column 
vector containing the expansion coefficients to be determined 
and v is the excitation vector, which corresponds to the tested 
incident fields. Next, we divide the computational subdomain 
into a number of contiguous sub-domains, and postulate that 
the solution on each subdomain can be found in the subpsace 
spanned by a number of precomputed subdomain solutions, 
referred to as the Macro Basis Functions and denoted by the 
column vector q;;,, Where 7 is the index of the subdomain and k 
in the index of an MBF defined on that subdomain. The above 
MBKFs need to be carefully chosen, and how to do this will be 
discussed in detail in Section IV. For the sake of simplicity, 
we will assume that the indices of all basis functions in a 
given subdomain are consecutive. It is then easy to identify 
blocks of the MoM impedance matrix associated with testing 
and basis functions on specific pairs of subdomains. One can 
also identify segments of the excitation vector, residual, or 
solution vector; those segments describe tested fields or current 
distributions on specific subdomains. If Q; denotes the matrix 
whose columns are comprised of the MBFs q,;,, then, for the 
ith segment of x, we assume that x; = Q;y,, for the i-th 
segment of x. The reduction of unknowns arises from the fact 
that vector y; contains much fewer unknowns than vector x; 
(typically by one to two orders of magnitude). By applying 
Galerkin testing, i.e., by choosing the set of macro testing 
functions identical to the set of macro basis functions, we 
obtain [6]: 

QF Z11Q; 


Q? Zin Qn yı Qty, 


QF Zy Qi Qu Zyn Qn Yn ARVN, 
where Q7 denotes the transposed conjugate of Q. In many 
implementations, just the transpose operation is applied, and 
it is difficult to say which one of these two options yields a 
better result. Since the matrices Q; have much fewer columns 
than lines, a very strong compression of the original system 
of equations is achieved. 

One should note that the original MBF approach [5] also 
employs elementary basis functions that “bridge” consecutive 
subdomains and that the real and imaginary parts of the MBFs 
are treated separately. These two refinements have either not 
been retained, or they have been integrated in different forms 
in subsequent MBF developments. 

It is interesting to note that the i-th block-line of (1) can be 
written as: 


QF (Z x]; — vi) =0 (2) 


MBF 


where x“™ is the solution obtained by using MBFs, and |g]; 
denotes the 7-th segment of a vector g (in the following, the 
brackets will in general be omitted). The expression between 
parentheses is nothing else than the opposite of segment 7 of 
the residual (r = v — Zx). This means that, as a result of 


Galerkin testing, the MBFs defined on a given subdomain are 
orthogonal to the segment of the residual corresponding to the 
same subdomain. 


II. HISTORICAL PERSPECTIVE 


As mentioned in the introduction, MBF-type methods 
rely on a divide-and-conquer approach to solve, through an 
integral-equation formulation, radiation or scattering problems 
involving structures that either have large electrical dimensions 
or fine features. The characterizing key features of these CEM 
frameworks are twofold: 

(i) Compression of the original MoM matrix equation by 
employing relatively few macro basis functions (MBFs) 
in order to exploit the low degrees-of-freedom (DoFs) 
that the physics-based equivalent current effectively at- 
tains, and reducing both the memory storage require- 
ments and solve-time significantly. 

Computation of the coupling between spatially (or spec- 
trally) distant MBFs in order to construct the reduced 
MoM matrix in a time-efficient manner. 


(ii) 


The objective of these CEM frameworks is to retain the 
low-order basis functions of high spatial resolution for the 
current (with minimum cell size /10) to be able to conform 
to arbitrarily shaped geometries, while reducing the DoFs for 
the current by employing MBFs. They present the additional 
advantage that existing MoM codes can be reused with only 
minor modifications. Within the MBF-type class of methods, 
one can recognize three widely-published CEM modeling 
frameworks. These are: 


e The Characteristic Basis Function Method 
(CBFM, [6]) which has been successfully applied 
to a large class of scattering [8], radiation [9], 
absorbing [10], as well as to waveguide and transmission 
line problems [11], [12]. Applications to planar antenna 
and microwave circuits have been described in [13]. 
This has been done typically by employing plane- 
wave-spectrum (PWS) generated CBFs for scattering 
problems (Sec. IV-A, and [14]), and primary, secondary 
or tertiary CBFs for radiation problems (Sec. IV- 
C), or a combination thereof [15]. CBFs partially 
overlap to preserve the continuity between electrically 
interconnected subdomains [16], and subdomain 
extension and windowing techniques are used to 
mitigate edge-truncation effects when generating CBFs 
on the interconnected subdomains (Sec. IV-A). CBF 
interactions for widely spaced subdomains have been 
computed rapidly, either using the Adaptive Cross 
Approximation (ACA) Algorithm (Sec. V-A, and [16]), 
an MBF-field interpolation technique (Sec. V-B, 
and [17]), or the Multilevel Fast Multipole Algorithm 
(MLFMA, [18]). The CBFM has shown to be highly 
paralellizable [19], [20], and a multilevel version of the 
CBFM has been described in [21], [22] and will be 
revisited in Sec. VII. 

The Synthetic-Functions Approach (SFX, [7]) applies 
the singular value decomposition (SVD) along with a 
thresholding procedure on the singular values to the 
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initially generated set of MBFs in order to orthonormalize 
and to retain only a minimum number of MBFs [23]. 
The SFX typically generates MBFs using point sources 
that surround the subdomain under excitation (Sec. IV- 
B). Furthermore, the SFX employs a separate and in- 
dependent set of low-order subsectional basis functions 
across the subdomain interfaces to electrically intercon- 
nect subdomains [24]. Far MBF interactions have been 
computed rapidly through an AIM approach (Sec. V-D, 
and [25]). It has primarily been applied to solve radiation 
problems [24], [26] and has also been hybridized with a 
multi-resolution approach [27]. 

The Macro Basis Function Method [5], which employs 
MBEFs obtained from both spectral (Sec. IV-D, and [28]) 
and spatial domain analyses. Both domains are also 
exploited for efficient computation of reaction integrals 
between distant MBFs (see for instance the multipole ap- 
proach both in Sec. V-C and [29], or the spectral domain 
approach in [30]). The method has been applied to both 
regular and irregular antenna arrays [31]. As mentioned 
in Sec. VI, closer inspection of iterative and MBF-based 
formulation has revealed an equivalence between specific 
types of MBF generation procedures and Krylov subspace 
iterative techniques, such as the Full Orthogonalization 
Method (FOM) [32]. Besides, a relationship has been 
established between the use of a block-diagonal pre- 
conditioner and the use of partially overlapping MBFs 
in iterative and iteration-free approaches [33]. 


Other more or less related subdomain-decomposition meth- 
ods are the Sub-Entire-Domain Basis Function Method 
(SED) [34], the Linear Embedding via Green’s Opera- 
tors (LEGO) technique combined with the eigencurrent ap- 
proach [35], [36], and a specific MBF domain decomposition 
technique, as described in [37]. 


IV. MBF GENERATION 


MBF-type approaches rely, for different subdomains, on 
an a priori choice of the subpsace in which the solution 
is expected to reside. This subspace is described by the 
MBEFs, whose choice is therefore crucial to obtaining accurate 
solutions. We describe below different methods that have been 
developed toward this end. 


A. Plane-wave spectrum 


The plane-wave spectrum approach [14], [38] calculates the 
current induced on each subdomain due to any electromagnetic 
field radiated by a source external to the domain. If the far-field 
condition is assumed, then the external field can be expanded 
in terms of a series of plane waves in the visible spectrum. 
Thus, according to the superposition principle, any induced 
current can be represented as a linear combination of the 
currents induced by the set of plane waves. For the m-th 
domain, the procedure can be mathematically expressed as 


die Za Pai (3) 


where P, is a matrix whose columns are the coefficients of 
the incident plane-wave field tested by the low-level basis 


functions in the m-th domain; Zmm is the impedance ma- 
trix comprised of the reaction terms between low-level basis 
functions in the m-th domain. 

In general, equation (3) is modified to consider extended 
subdomains [14], [38]. The purpose of this extension is 
twofold. First, the edge effect due to the domain truncation 
is moved away and, second, it enables us to include the near- 
field contributions of the region closest to the domain. 





Fig. 1. 


MBFs generation based on the plane-wave expansion. 


Fig. 1 illustrates this approach. The continuous thick trace 
shows the boundary of the domain in which the MBFs are 
being generated whereas the dotted thick trace delimits the 
extended domain wherein the currents induced by each plane 
wave are calculated. 

After discarding the currents in the extension, the computed 
currents are filtered using the singular value decomposition, 
which yields the final set of MBFs and guarantees the orthogo- 
nality between the MBFs. Thus, the SVD entails the following 
matrix factorization: 


Jin = O,=.v; (4) 


where the diagonal of the matrix X, contains the singular 
value of the decomposition. The final MBFs coefficients Qm 
are calculated by retaining only the columns in Q,, whose 
normalized singular value is above a prescribed threshold 7. 
Hence, the i-th column is only retained if o;/0, > T. Typical 
values for this threshold ranges from 107? to 1075. 

The plane-wave spectrum approach typically yields a higher 
number of MBFs as compared to the previous approaches. 
However, the computed MBFs do not depend on the excitation. 
Consequently, it is usually preferable to analyze problems 
that involve multiple excitation sources, as for example in 
monostatic RCS computations. A modification of this approach 
is to employ spherical waves instead of plane waves, as 
suggested in [10]. 


B. Point sources 


Another approach for generating the MBFs consists of 
replacing plane waves by a number of point sources distributed 
over a given surface that surrounds the subdomain of interest 
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[39]. This approach may be regarded as relying on the surface 
equivalence principle, according to which the sources external 
to the surface can be replaced by equivalent electric and 
magnetic currents on the surface [40]. In addition, for spherical 
surfaces, the equivalent current can be limited to electrical 
currents only. Such equivalent currents, in principle, ensure the 
completeness of the MBF basis formed in this way. The main 
reason why the base may not be truly complete in practice 
is the limited sampling of the equivalence surface. This may 
become an issue when the surface very closely wraps the 
subdomain of interest, and the approach becomes virtually 
impractical when the subdomains are connected, since it then 
becomes difficult to let the equivalence surface partition the 
subdomains, unless the equivalence surface entirely includes 
the extended subdomain introduced in the previous section. 


C. Primary and secondary MBFs 


Another way of generating the MBFs uses primary and 
secondary current distributions [9]; it is particularly suitable 
(but not limited) to the analysis of antenna arrays. Indeed, 
for the analysis of mutual coupling in arrays, it is generally 
sufficient to provide all the embedded element patterns as 
well as the array impedance matrix. The above quantities can 
be obtained from the solutions derived by exciting the array 
at each of the individual ports. The construction of MBFs 
may then be obtained from the excitation of the antenna in 
isolation, followed by the excitation of the other elements 
by the fields radiated by the first element. From the MoM 
point of view, this solution is obtained by multiplying blocks 
of the MoM impedance matrix. More precisely, the primary 
MBF on domain i corresponds to fp; = Ze v;i, where v; 
is the excitation vector on antenna (or subdomain) 2, while a 
secondary MBF on antenna 7 is obtained from the equation 
fsa = Zi; Zjifp i- In order to enrich the set of MBFs, it 
is logical to employ the primary and secondary MBFs on all 
antennas, by considering every possible excitation, or at least 
secondary MBFs created from primaries on every neighboring 
subdomain. This approach usually provides excellent results 
on arrays of disconnected elements. For arrays of connected 
elements, combining this idea with the use of extended sub- 
domains, has proven to be very efficient and accurate, as has 
been explained in Sec. IV.A. 

The idea of primary and secondary MBFs can be extended 
to higher multiple-scattering orders, by including the tertiary 
MBFs as done for instance in [41], [42] and [43]. There is 
virtually no limit to the orders that can be considered, albeit 
at an increased computational cost. As explained in Section VI, 
the completeness of MBFs bases can, in principle, be achieved 
by considering virtually unlimited orders (i.e., only limited by 
the number of unknowns in the problem), though this is not 
a viable option in practice. Fortunately, very high accuracy 
can be achieved with orders limited to 2 or 3, especially when 
extended subdomains are used. In [33], the somewhat complex 
process of domain extension has been made implicit by first 
pre-conditioning the system of equations. The preconditioner 
utilized is a nearest-interaction preconditioner, which can be 
regarded as an extension of the shield-block preconditioner 


introduced in [44]. In a nutshell, an extension Sja is associated 
with each subdomain S; in this approach (see Fig. 1). In the 
following, segments of vectors and blocks of matrices will 
be associated with different subdomains and their extensions 
by using indices 7 and 7°, respectively. The preconditioned 
system of equations then reads Z x = w, with the following 
definitions: 


with v; and v;a corresponding to segments of the original 
excitation vector, and 


P; = Zijie Zija (6) 
Y; = (Zin—P; Zia i) t (7) 


There are two reasons for doing this. First, the preconditioned 
system of equations ensures faster convergence of Krylov- 
based iterative techniques. Second, MBFs of order n can be 
obtained simply through multiplication to the left of a primary 
MBF by a number of consecutive matrices. By “consecutive” 
we mean that a matrix with first index k must be multiplied 
to the left by a matrix with second index k. In Section VI, it 
will be shown how such MBFs can be combined to construct 
Krylov subspaces, in which solutions are sought in iterative 
schemes. 


D. ASM-MBF 


The ASM-MBEF approach is limited to regular arrays of 
antennas or scatterers [45]. For array problems, one seeks the 
solutions (current or field distributions) over the entire array 
when an arbitrary element is excited. Therefore, it makes sense 
to obtain the MBFs from the field or current distribution in an 
infinite array when a single element is excited. As explained 
in [46], this problem can be solved as the superposition of 
infinite-array problems (with all elements excited) by scan- 
ning through every possible inter-element phase shift. In one 
dimension, this is expressed as: 


2T 
in= a | TOW ay O) 
2m Jo 
where Jan is the current on element m when element O is 
excited, and J® (ọ) is the infinite-array current at the same 
position within the unit cell, for an inter-element phase shift 
equal to %. By superposition, the current distributions obtained 
on successive elements when a single element is excited forms 
an excellent basis for an arbitrary excitation law. Even the 
effects of array truncation can be well represented in this 
basis, since currents “reflected” by the edges of the array 
may form current distributions that are very similar to those 
obtained from “direct” waves launched by a single element 
in the infinite array. This might not hold true for elements 
located right at the edges (or corners) of the array, in particular 
when the elements are complex and connected with each other. 
Therefore, to improve the accuracy, a few current distributions, 
obtained in a 2x2 array, are added to the set of MBFs. 
It has been found that this approach leads to a very fast 
convergence of the solution w.r.t. the number of points used 
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to discretize the ASM integration, and excellent accuracy has 
been realized using about 20 MBFs per element. More impor- 
tantly, for the reasons explained above, the MBFs obtained 
in this manner are excitation-independent. An open-source 
Matlab software for the example of linear dipole arrays has 
been described in [47]. An extension of this methodology for 
arrays of plasmonic rods has been detailed in [48]. 


V. FAST MBF INTERACTIONS 


The construction of the reduced matrix equation in (1) 
requires us to compute many blocks of the form 


A 
Qr Zinn Qn 


and it is desirable to perform this computation in a time- 
efficient manner. From a physics point-of-view, the factor 
Zinn Qn represents the excitation matrix Vmn due to the 
source MBFs on the nth subdomain, whose radiated E-fields 
are tested on the mth subdomain. As the source and obser- 
vation subdomains become electrically well-separated in free 
space, the DoFs of any such subdomain excitation vector (col- 
umn of Vmn) reduces. In fact, for extremely large separation 
distances, each excitation vector practically represents a single 
incident plane wave field (thus only one DOF, or one mode). 
One can exploit this phenomenon to rapidly compute (10), 
either through a field expansion method employing only the 
first few dominant modes, or by using an algebraic method 
exploiting the low-rank nature of Zmn. 
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Fig. 2. Normalized singular value spectrum logy 9(|on|/|o1|) of the coupling 
matrix block Zmn between a pair of 2A x 2A plates (A/10 meshing, 2320 
RWGs in total), as a function of their separation distance d when (a) facing 
each other, and; (b) in a side-by-side configuration. 








As regards the rank-deficiency of Zmn, Fig. 2(a) and (b) 
show how the singular value spectrum of Zmn depends on 
the separation distance d between a pair of 2\ x 2A plates. 
Note that, when defining the effective numerical rank as r = 
rank(Zmn) = |o,|/|o1| = 107?, i.e., when r is the number of 
singular values that are within 107? from the largest singular 


value, one can observe how r decreases as a function of d 
[cf yellowish region in Fig. 2(a)]. For d = 10A, the effective 
rank is r ~ 10, which is less than 0.5% of an equally large 
full-rank matrix. It is also observed that the effective rank 
decreases even faster for plates that are placed side-by-side. 
For instance, for d = 10A, we find that r ~ 4, which is smaller 
than 0.1% of an equally large full-rank matrix. Clearly, both 
the subdomain sizes and orientation play an important role in 
the degree of rank-deficiency of Zmn and, consequently, on 
the computation time of the reduced matrix elements in (1). 


A. The Adaptive Cross Approximation (ACA) algorithm 


The Adaptive Cross Approximation (ACA) algorithm, origi- 
nally developed by Bebendorf [49], approximates the Nm x Nn 
rank-deficient matrix block Zmn through the low-rank block 
factorized matrix Zmn = U,,%™*'V,,"*%". This is advan- 
tageous because (10) can then be computed rapidly using a 
minimum number of multiplications as 


(Q7U,,)(V,, Qn). (1) 


A very important feature of the ACA algorithm is that the 
matrices Um and V,, are constructed on-the-fly, without a 
priori knowledge of the entire original matrix block Zmn; 
the iterative ACA algorithm dynamically selects certain rows 
and columns of Zmn and, in conjunction with a normaliza- 
tion procedure, these normalized vectors form the rows and 
columns of the matrices V,, and Um, respectively. Indeed, for 
well-separated groups of RWGs (Rao-Wilton-Glisson double- 
triangle basis functions), the electric field at the observation 
group m produced by any source RWG can be expressed 
as a linear combination of the fields produced by only a 
few of these source RWGs (source sampling). Likewise, the 
electric field tested at the observation group m produced by 
any source RWG can be expressed as a linear combination 
of the fields tested by only a few of these observation RWGs 
(field sampling). Hence, a cross-approximation technique can 
be used to adaptively construct the subsets of relevant source 
and observation RWGs. 

The ACA algorithm is purely algebraic in nature, easy 
to implement, and can be used irrespectively of the kernel 
of the integral equation, basis functions or type of integral 
equation formulation. The ACA algorithm has not only been 
applied to solve low-frequency EMC problems [50], but also to 
solve electrodynamic antenna problems involving oscillatory 
kernels using an MBF approach [16]. Since the ACA algo- 
rithm approximates Zmn through the product UmVn, most 
of the non-selected elements of Zmn are predicted through 
linear interpolation, i.e., from the product U,,,V,,; hence, the 
time-harmonic nature of the fields is not accounted for. The 
ACA algorithm may therefore require more iterations than 
a more physics-based approximation technique, such as the 
multipole approach as explained below. Also, the compu- 
tational overhead of the ACA algorithm becomes excessive 
for a relatively large numerical rank of Zmn (e.g., for small 
subdomain separation distances), so that a direct element-by- 
element computation of Zmn is more efficient. The interested 
reader may refer to [50], where a pseudo-code of the ACA 
algorithm in Matlab notation can be found. 
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Fig. 3 shows the matrix fill-time of the ACA algorithm 
for building the matrices Um and V,, — when applied to the 
case shown in Fig. 2(a) — relative to the time needed when a 
direct matrix filling approach is used for building Zn. It is 
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Fig. 3. Fill time of the ACA constructed matrix blocks Vn and Um, relative 
to a full element-by-element filling approach of Zmn, as a function of the 
separation distance d for a pair of plates facing each other [cf. Fig. 2(a)]. 


evident from Fig. 3 that the speed advantage of the ACA over 
a direct matrix filling technique is significant. For instance, 
for d > 0.5A, the ACA algorithm requires less than 15% of 
the time needed to fill a full MoM block on an element-by- 
element basis. This is true even for ACA thresholds as low as 
1078, which means that the relative ACA approximation error 
[Vn Um — Zmalle/||Zmnllp < 1073, where ||- ||, denotes the 
Frobenius norm. Hence, for electrically large problems, the 
average ACA matrix fill time typically takes only a few percent 
of that needed in a direct matrix filling approach. As an 
alternative to the ACA technique, matrix compression based 
on the incomplete QR decomposition [51] has been used in 
[33]. 


B. Tested field interpolation 


Another technique, which also exploits the DoF of the 
field radiated by the MBFs, is based on the conventional 
interpolation of the radiated field [17]. Thus, the tested field 
can be computed by calculating the field in a small grid over 
the observation domain and, then, retrieving the field in the 
low-level basis functions via interpolation. Fig. 4 shows this 
interpolation scheme for a planar geometry. 

In order to rapidly compute the field radiated by each MBF 
in the n-th source domain over the m-th observation domain, 
the matrix Vmn relating the coefficients of the source domain 
MBFs and the p-component of the field (p = x, y or z) in the 
interpolation grid are computed. By invoking the reciprocity 
theorem, the entries of this matrix can be expressed as [17]: 


Vnnli.J] = J fng Ex(tma)dS, 


wherein Tm, is the i-th observation point in the interpolation 
grid for the m-th MBF domain; E’7(r) is the field radiated 


(12) 


Interpolation grid Extended observation 








Fig. 4. 


Interpolation grid for fast computation of the reduced matrix. 


by an infinitesimal dipole at the spatial point r and oriented 
along the p-axis; and, f„ j is the j-th low-level basis function 
in the n-th MBF domain. Once these matrices have been 
computed, the field in the sampling grid can be calculated by 
post-multiplying with the MBF coefficients Egrid =VinnQn. 
Thus, the need for computationally expensive integrations in 
the source domain is obviated. It is remarkable to note that 
this approach is compatible with more advanced interpolation 
schemes such as those proposed in [52] wherein a phase 
extraction is carried out first to further reduce the DoFs. 

The interpolation scheme is illustrated by means of the 
example shown in Fig. 4 where the bistatic analysis of two 
square plates with edge lengths of 2A is considered. Both 
plates lie in the same plane and the distance between them is 
1.5. The frequency is chosen to be 300MHz. For this case, 
the interaction between both blocks is calculated by using a 
9x9 interpolation grid, which enables one to reduce the time 
to compute the reaction term between the MBFs of both plates 
from 4.28s to 0.37s. We note an excellent agreement in the 
entire dynamic range of the bistatic RCS. 











Conventional CBFM 
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Fig. 5. Bistatic RCS between two square plates with a distance of 1.5 and 
edges equal to 2A. 


Another approach involving interpolation for estimating 
MBF interactions is described in [53]. It has been developed 
for the analysis of irregular arrays of identical antennas (or 
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scatterers) and it produces a very simple model for the 
interactions between MBFs versus relative position. It is based 
on three physical transformations: far-field extraction, phase 
removal and change of the distance variable. In this way, a 
harmonic-polynomial model is obtained, valid for any relative 
position in the plane, through the explicit calculation of 
interactions at a few tens of relative positions. 


C. Spectral approaches 


Interactions between subdomains can be speeded up by ex- 
ploiting integral representations of the scalar Green’s function. 
This is particularly fast when, in such representations, the 
dependence on source and observation coordinates is sepa- 
rable. In practice, the separable form is generally identified 
with a plane wave, expressed by a complex exponential. 
Two categories of spectral approaches have been developed 
in the literature. The first one is associated with multipole 
decompositions, while the second one is associated with waves 
radiated from a given reference plane. 

a) Multipole approach: The derivation of this approach 
is provided in [29]. The final result reads as follows. If Ë, is 
the radiation pattern of a conjugated Macro Testing Function 
and F, is the radiation pattern of a conjugated Macro Basis 
Function, the interaction between them can be written: 


ref] F* . F, T(k,?, ù) dU 


where T (å) is the translation function appearing in multipole 
decompositions, within a constant factor, k is the free-space 
wavenumber and 7’ is the vector distance between reference 
points of the source and observation domains. The integration 
domain U corresponds to the unit sphere, to which the unit 
vector ù points. This approach allows the computation of 
the interactions between subdomains without computing the 
interaction matrix Z;;. The only constraint is that the distance 
between subdomains should exceed a certain minimum, whose 
value is of the order of half wavelength. Fig. 7 illustrates 
the accuracy of the multipole-based method with 40x40 
integration points over the unit sphere, for the antenna shown 
in Fig. 6. The MBF considered is a primary (direct excitation 
of one antenna); the solid line provides the magnitude of the 
interaction versus distance in wavelengths, while the dashed 
line povides, on the same log-scale, the magnitude of the 
difference between results obtained using the MoM matrix 
approach described in Section 2 on one side, and the multipole 
approach on the other side. It can be seen that the quality 
suddenly degrades for very small distances. However, this 
sudden change happens when the antennas are nearly touching 
each other. More precisely, if the acceptable threshold is 
defined at a 1 % error level, then for the 5 cm wavelength, the 
tip-to-tip distance between antennas should be at least 0.5 cm, 
while that distance is only 0.2 cm for the 2.5 cm wavelength 
case. 

If N is the number of elementary basis functions on a given 
subdomain, the complexity of computing interactions between 
the MBFs is typically reduced from N? to N. Assuming 
sudomains of the order of one wavelength and a relatively 
coarse mesh; the time saving is smaller for larger domains 


(13) 


and larger for finer meshes. Such time-saving has been demon- 
strated in the case of arrays of broadband conducting antennas 
in [29], and has been extended in [32] to subdomains made 
of penetrable bodies. An extension to printed antennas is 
described in [54]. In the latter case, the Green’s function 
is decomposed into a spherical wave, related to an average- 
medium term, as well as into cylindrical waves [55]. In both 
cases, the MBF interactions are computed using multipoles. 
The treatment of the terms related to the cylindrical waves 
has been described above. For the multipole-based treatment 
of the terms related to the cylindrical waves, the complexity 
is proportional to N!/? where N is the number of elementary 
basis functions per antenna, and to the number of cylindrical 
waves needed, which is typically in the order of 10. For 
this case, the computation time for printed structures is only 
marginally larger than what it is when the subdomains are 
interacting in free space. 
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Fig. 6. Discretization of the bowtie antenna considered in the multipole 
analysis of figure 7. 
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Fig. 7. Interactions (solid) between primary MBFs defined on a pair of 
antennas versus center-to-center distance (vertical shift in Fig. 6). Dashed: 
error incurred by multipole-based approach. Top: 5 cm wavelength. Bottom: 
2.5 cm wavelength. 


b) Waves from a reference plane: Assuming a reference 
plane XY, the scalar Green’s function can be written as a 
continuous spectrum of plane waves, characterized by their 
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lateral wavenumbers ky and ky. If k2 + k? < k?, then the 
plane wave is propagating along z, otherwise it is evanescent. 
Using such a decomposition, the interactions between MBFs 
can be written as 


r=x |J fy eI (Fe Aethy Av) £* Gilkey, ky) dka dky (14) 


where f; and fẹ are the Fourier transforms, or “patterns”, of 
MTFs and MBFs in direction (kz, ky, kz)/k, which becomes 
complex outside the unit circle in the (kz, ky)/k plane; Az and 
A, are distances between reference points of the subdomains 
in the XY plane; G is the spectral-domain representation 
of the dyadic Green’s function and K is a constant. Such 
an approach, specialized to analytically-derived CBFs, has 
been presented in [56]. Since the MBFs are defined over 
domains that are substantially larger than those of elementary 
basis functions, their pattern is relatively narrow; hence the 
integration domain in wavenumber space can be strongly 
reduced. Examples of this approach are given in [30] in the 
case of printed antennas. 


D. FFT-based approach 


In [57], the interactions between MBFs are obtained using 
the AIM for printed structures. This may be viewed as a 
fast spectral approach, since the space-domain convolution 
between MBFs, MTFs and the Green’s function are written 
as products in spectral domain. In that approach, forward and 
backward 3D FFTs are exploited to compute the space-to- 
spectral and spectral-to-space domain transforms. This may 
be regarded as one of the most effective MBF-interaction 
approaches to date. Reference [57] also provides expressions 
for the complexities of the different interactions techniques 
as well as validations for large problems, such as arrays of 
printed antennas. 


VI. RELATION WITH KRYLOV SUB-SPACE ITERATION 
A. Reformulation of Krylov iterations 


As reminded in Section II, the concept of Macro Basis Func- 
tions can be expressed in a relatively compact form. Krylov 
subspace techniques, essentially developed in the seventies, 
are also based on a few key ideas, some of which also appear 
in MBF methods. 

However, the use of specific tools, as for instance the 
properties of Hessenberg matrices resulting from the Arnoldi 
orthogonalization process [58], may somewhat obscure the 
basic ideas behind Krylov methods. In this section, we propose 
an alternative —or perhaps simplistic— formulation for two 
popular Krylov-based methods, namely the FOM (Full Or- 
thogonalization Method) and GMRES (Generalized Minimal 
Residual), both of which are described in the seminal paper by 
Saad and Schulz [59], published in 1986. As explained further 
below, that alternative formulation may incur a very marginal 
reduction of efficiency. However, the formulation proposed 
here should further clarify the relationship between Krylov 
subspace iterative techniques and the MBF approach. This will 
be explained below in two steps. 


First, let us assume a relatively well preconditioned system 
of equations A x = b, with A sufficiently close to a unit matrix. 
The reader is referred to [58] for more precise figures of merit 
of preconditioning and to [60] for recent advances regarding 
preconditioning in the framework of integral-equation solution 
in high-frequency electromagnetics. For the above system of 
equations, with xg as an initial guess, the first residual is 
ro = b — Axo and the simplest possible iteration [58] is 
obtained by considering at iteration k a correction equal to the 
residual rg—1; hence x, = Xķ—ı +b — Ax,_;. As compared 
to Xk—1, Xz has high chances to be closer to the exact solution 
Xk—1 + AT! (b — Axķk—1ı), because A is relatively close to a 
unit matrix, as a result of preconditioning. The convergence 
of this procedure is dictated by the eigenvalues of the iteration 
matrix (I — A), which are unfortunately not known a priori. 
As has been very well summarized in [61], Krylov iteration 
essentially consists of keeping all the approximants obtained 
up to the step k and to recombine them to obtain a more 
accurate solution. This is equivalent to searching x’ = x — xo 
in the subspace spanned by the successive residuals ro, r1, etc. 
Based on the above, it is easy to prove that this subspace can 
also be written as Span{ ro, Aro, A” ro, eh ro}, which 
is the Krylov subspace of order k. Each of the vectors 
sld) = AT ' ro describing this subspace, which we will name 
the generating vectors, may be regarded as MBFs spanning 
the whole computational domain. In the following, we will 
denote by Q the matrix whose columns are formed by using 
the consecutive generating vectors. 

Second, one needs to establish a set of conditions which will 
determine the scalar coefficients yo...y,—1 that multiply each 
generating vector in the final estimate of x, so that we can write 
x’ = x— Xo ~ Qy. In passing, to avoid dealing with the initial 
guess Xo, the system of equations may be rewritten as Ax’ = 
ro. A simple approach for finding the vector of coefficients y 
consists of testing the original system of equations with the 
generating vectors, i.e by using 

Q” AQy =Q" ro (15) 
This is mathematically equivalent to the Full Orthogonal- 
ization method (FOM, [59]), which imposes the residual to 
be orthogonal to the Krylov subpsace, and hence to each 
generating vector s“ composing Q. This also has the same 
form as an MBF approach performed over a single domain. 
The reduced system of equations obtained in this way will in 
general be ill-conditioned, because of quasi-linear dependence 
between generating vectors. This is why it is necessary to 
orthogonalize the generating vectors composing Q. This can 
be achieved using the Arnoldi procedure, inspired from the 
Gram-Schmidt method. As compared to (15), an alternative 
approach consists of testing the initial system of equations 
with the generating vectors, each left-multiplied by matrix A. 
This may be written as 


(AQ)” (AQ)y = (A Q)” ro (16) 


which may be viewed as the normal equation which minimizes 
the residual Ax’ — ro in the least-squares sense. This is 
equivalent to the GMRES approach [59]. 
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In practice, the above methods can be implemented as 
follows. As the consecutive generating vectors A” ro are 
created, they are orthogonalized and placed into matrix Q. 
At every step, the following system of equations is solved: 
B” AQy = BË ro, with B = Q for FOM and B = AQ 
for GMRES. The solution is then given by x = xo + Qy, 
The A Q vectors can be computed at very low cost, based on 
operations already carried out to build the Krylov subspace. 
This observation makes it possible to use only one matrix- 
vector product per iteration instead of two. The simplicity of 
this new formulation may incur a slight extra computational 
cost. Indeed, standard implementations of FOM and GMRES 
make use of the Hessenberg matrix [59] resulting from the 
Arnoldi procedure used to orthogonalize the Krylov subspace. 
It is an upper-triangular matrix with an extra non-zero diagonal 
below the main diagonal (it is related to the R matrix resulting 
from the modified Gram-Schmidt procedure). Consecutive 
solutions of the larger problems obtained at each iteration 
can be accelerated by keeping in memory previous pivoting 
operations on the growing Hessenberg matrix. However, as 
long as the number of iterations is much smaller than the 
number of unknowns (say by at least one order of magnitude), 
the solution time of the reduced system of equations is not 
the limiting factor. The dominant part of the computational 
effort is associated with the product between matrix A and 
consecutive generating vectors. Under those circumstances, 
the simple procedure proposed here is, in practice, equally 
efficient. 


When the number of iterations becomes large, the gener- 
ating vectors may require too much memory and their or- 
thogonalization may become too expensive (linear growth per 
iteration). Also, they may lead to poor conditioning because 
they are only marginally linearly independent, a problem that 
can get exacerbated by numerical roundoff. Then, the iteration 
may be restarted, considering the solution obtained after a 
number of iterations as the new first guess. 


B. Relation between Krylov iteration and MBF approach 


Regarding MBF construction, a particularly well-posed ap- 
proach consists of using ”primary” and ”secondary” MBFs 
[9], extended to higher orders in [41]. In a nutshell, the current 
distribution obtained on a given subdomain impresses fields on 
another subdomain, in which currents are induced when that 
other subdomain is taken in isolation; this process is continued 
in a multiple-scattering approach. Let us denote by S; the 
subdomain of interest and by S; all the other subdomains. If 
the currents on other subdomains are known exactly, then the 
MBFs they induce on subdomain S; form a complete set, with 
known coefficients. This can be proven from the t-th block- 
line of the system of equations, which represents the testing 
of fields on S;. If the Z;; matrix denotes blocks of the system 
matrix, x; denotes solutions on Sj and v; the segment of 
the excitation vector standing for testing of incident fields on 
Si, then the multiple-scattering process produces the following 


general solution on Sj: 


X; = Ze: Vi -X 052i Xj (17) 


j+i 


where the Z;;'Z,, x; terms are the secondary MBFs and the 
a; coefficients are yet to be determined. It is also obvious 
from the 7-th block-line of the system of equations that an 
exact solution for x; can be obtained with a; = 1 for all 7’s if 
the x;’s are known a priori. Of course, this condition sounds 
difficult to satisfy; however, in common with the Krylov-based 
approaches, the fact of keeping free all the coefficients that 
multiply the generated MBFs provides important degrees of 
freedom (DoFs). To a large extent, those DoFs may compen- 
sate for the deficiency of working with MBFs generated from 
inaccurate current distributions x; on the subdomains S;. More 
DoFs are obtained by adding higher-order multiple-scattering 
MBFs, as explained in Sec. IV.C. 

Two challenges appear when implementing the above pro- 
cedure. First, in a multiple-scattering process, the number of 
generated MBFs increases exponentially. Second, connected 
subdomains may lead to non-physical MBFs on S;, with 
nearly-singular current distributions along the contour of S;. 
To circumvent the latter problem, several authors [62], [63] 
opted to extend subdomain S; (see also Sec. IV.A), with 
a connected auxiliary subdomain S and to retain as an 
MBF the current induced only on subdomain S; (see Figs. 
1 and 2). As already mentioned in Section IV, in [33], it is 
proven that this approach is equivalent to the classical MBF 
approach (i.e., without subdomain extensions), provided that 
the system of equations is modified using a nearest-interactions 
preconditioner, which may be viewed as an extension of the 
shielded-block preconditionner proposed in [44] for discon- 
nected periodic structures. We will denote by Zx = w the 
system of equations preconditioned in this manner. Details 
about the combination of this preconditioning with compres- 
sion techniques (based on ACA or QR) have been provided in 
[48]. Given such preconditioning, the extension of subdomains 
becomes implicit, and multiple-scattering MBFs are simply 
generated via multiplication to the left by consecutive blocks 
of the preconditioned system of equations. In [33], it is proven 
that, if for a given excitation all multiple-scattering MBFs are 
generated up to order p, then the Krylov subspace of order 
q < pcan be constructed with the exclusive help of the MBFs. 
Mathematically, this reads 


a =O thy (18) 


where s is segment i (corresponding to currents in the 
subdomain 7) of generating vector of order q < p (i.e., the q-th 
vector defining the Krylov subpsace, Sec. VIA); the columns 
of QP ) are the MBFs on S; up to order p; and fp, contains 
coefficients that can be unambiguously determined. Besides 
the mathematical proof given in [33], an intuitive argument 
may be advanced to justify the above statement, namely if the 
generating vectors of the Krylov subspace are created through 
consecutive left-multiplications by the entire matrix Z, then 
they can only involve series of “consecutive” blocks of Z, 
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left-multiplying a segment w; of the excitation vector. As a 
reminder “consecutive” means here that, if block Z;; left- 
multiplies block Z;;, then j = k. If the preconditioned system 
of equations is used, the segment w; actually corresponds 
to a primary MBF, while the products that appear in the 
formation of the generating vectors simply correspond to 
multiple-scattering MBFs. Hence, the segments of the different 
generating vectors are linear combinations of MBFs, and as 
far as the Krylov subspace is considered as a meaningful 
basis for the entire solution, the above property indicates that 
the multiple-scattering MBFs are going to form an equally 
meaningful basis, when used on their respective subdomains. 
This alleviates the problem with the issue of completeness of 
the MBF subspace. However, as mentioned above, the number 
of MBFs generated in this way rapidly becomes prohibitive, so 
much so that higher-order MBFs or MBFs generated through 
interaction between very distant subdomains in general need 
to be discarded. If deemed necessary, they can be replaced by 
MBFs generated by using a different approach. For instance, 
it is clear that MBFs generated via distant interactions can be 
very well represented by MBFs generated using plane waves, 
such that MBFs of orders larger than two may be created based 
only on interactions between contiguous subdomains. 
Regarding the condition imposed by the MBF approach 
to obtain the solution of Z x w, another link can be 
established with the FOM Krylov iteration [33], which we 
have re-formulated in Sec. VIA. Using the definitions in 
Section II, we provide here a derivation of that property that is 
more direct in comparison to that given in [33]. Let us denote 
by rP ) the i-th segment of the residual vector, obtained from 
the multiple-scattering approach up to the order p, while all 
MBFs are retained. By construction, the MBF solution is such 
that rP ) is orthogonal to the MBFs defined on S; (see (2) and 
comment below), i.e., QP” rP) = 0. Hence, using (18), we 
have, for q < p, sD H pP) a QPF rP) = 0. The same 
reasoning can be held for all subdomains S;, such that: 


y sD H rP) = sD H r) — 9 


(19) 


because each term of the sum is zero, which proves that 
the entire generating vector s‘% is orthogonal to the entire 
residual rP), In other words, the MBF approach satisfies the 
orthogonality conditions that characterize the FOM solution of 
order q equal to or smaller than the multiple-scattering process 
p (for the FOM condition, see (15) and its interpretation 
below.) This however supposes that all multiple-scattering 
MBFs are kept up to order p, which is not truly practical 
(see above). Also, the MBF solution of order p satisfies more 
conditions, at the cost of having to work with a larger system 
of equations, with more unknowns to be determined than is 
needed in the FOM approach of same order. 

Now let us examine an important question regarding the 
comparison of performance, in terms of accuracy and compu- 
tational cost, between Krylov-based iteration and MBF-based 
solution. In [33], an empirical comparison has been performed 
on different types of examples (arrays, spheres and aircraft). 
This comparison is limited to the following conditions: the 
MBFs are generated in an excitation-specific way and through 
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a truncated multiple-scattering process, only up to order two 
(with the nearest-neighbor preconditioner the higher orders do 
not seem to significantly improve the accuracy). Numerical 
experiments have led to the following (albeit preliminary) 
conclusion: for equal cost in terms of computation time and 
memory, the accuracy of both methods is similar when the 
number of iterations of the FOM approach is equal to the 
number of MBFs per subdomain. This just appears to be a rule 
of thumb and needs to be further tested with other examples. 
We provide below one more example, involving a connected 
array of bowtie antennas, also studied in [64]. The 5x5 array 
is shown in Fig. 8. In the preconditioning step, all of the 
neighboring elements are included in the auxiliary subdomain. 
Fig. 9 shows the port currents obtained when only the element 
1 is excited (top line). The ports do not contain a series 
impedance; simulations are carried out for a frequency of 10 
GHz. It can be seen that the port currents on the other elements 
are not very low as compared to the one in the excited element, 
as a result of extremely strong coupling. Thus, this strongly 
coupled array forms a good test case for numerical methods. 
On the same plot in log scale, the differences between exact 
and approximate solutions are shown for the MBF approach, 
with 9 MBFs per subdomain (corresponding here to one 
antenna) for the preconditoned system of equations, and the 
errors obtained for the FOM approach with also 9 iterations. 
It can be seen that a comparable error level is achieved, with a 
slight advantage for the MBF approach. Quasi-identical error 
levels, within 1.7 dB on the average, were obtained with 11 
iterations for FOM, instead of 9. This new example supports 
the rule of thumb referred to above. Further tests are now 
being carried out in the field of metamaterials. 

Given this rule of thumb, one may wonder about the real 
advantage of using MBFs, as compared to iterative techniques, 
and we offer at least four. First, it is obviously advantageous 
to reduce the number of DoFs when, because of the level 
of geometrical detail, the discretization of the structure needs 
to be much finer than the usual \/10 condition. This is 
particularly true for antenna applications, where the often 
complex feeding region needs to be modeled with many 
elementary basis functions. Second, when appropriate general- 
purpose MBFs can be found, the efficient solution for multiple 
right-hand sides offers an important advantage. An example 
where such MBFs can be easily found is non-periodic arrays 
of antennas; though more challenging structures involving con- 
nected elements may also belong to that category. Third, being 
non-iterative, the MBF approach is much easier to parallelize. 
Finally, we also observed for the case of scattering by spheres 
that the rule of thumb referred to above tends to break down 
for resonant structures, with an important advantage, in terms 
of accuracy versus computational resources, in favor of the 
MBE approach. This seems also to be the case for extremely 
finely meshed structures, such as those studied in [48]. 

A further link between Krylov-based methods and the MBF 
approach may lie in the “restart” procedure [59]. For Krylov- 
based iterations, this means that the solution at a given point 
can be considered as the new first guess for the iterative 
technique. Similarly, the MBF-based solution obtained on a 
given subdomain can be transformed into just one MBF, while 
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Fig. 8. Mesh used for simulation of 5x5 array of connected bowtie antennas. 
Element spacing: 1.044 cm horizontally and 1.00 cm vertically. 
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Fig. 9. Exact solution for port currents (numbered left to right, top to bottom, 
see Fig. 8), and errors obtained with 9 MBFs per subdomain and with 9 FOM 
iterations. 


new MBFs are added, based either on plane-wave excitations 
or on higher-order multiple scattering. This enables one to 
refine the solution without augmenting the dimension of the 
final system of equations. A preliminary study of this approach 
has been carried out in [64] and may lead to a new direction 
of research on MBFs and CBFM. 


VII. MULTI-SCALE MBF ANALYSIS 
A. Degrees of freedom 


Methods relying on MBFs benefit from a reduction in the 
number of unknowns. This reduction offers a large number 
of advantages such as the well-known memory saving or the 
use of conventional MoM techniques that have been typically 
limited in their application to electrically small problems 
(e.g., [65]). As was previously detailed in Section H, MBFs 
are defined over a set of contiguous low-level basis functions. 
In the analysis of antenna arrays, it becomes natural to apply a 
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partitioning so that each subdomain corresponds to an element 
of the array (e.g., [6]). However, there is nothing that prevents 
us from defining MBFs for a group of adjacent elements. For 
the case of electrically large antennas or scatterers, there is no 
straightforward strategy for the partitioning and the preferred 
choice is to select the subdomains by grouping basis functions 
inside certain canonical geometries (e.g., cubes) [38], [66]. 

At this point, one might wonder how large the size of 
the subdomains should be to define the MBFs. From the 
solve-time point of view, some criteria have been recently 
proposed [67] based on minimizing the complexity. However, 
it does not consider the reduction of the number of unknowns, 
which is one of the key features of the MBF approach. 
Although the number of degrees of freedom (DoFs) of radiated 
and scattered fields is well-known [68], [69], to the best of 
authors’ knowledge, that is not the case when considering the 
number of DoFs (and so the number of MBFs) for the currents 
induced on an arbitrary surface. 

To further examine this issue, connected with multilevel 
MBE approaches, we will now present some numerical simu- 
lations. Let us first consider a given geometry discretized by 
means of low-level basis functions, i.e., RWG functions. The 
object is illuminated by a set of plane waves. As described in 
Section IV, the current induced by any arbitrary incident field 
can be calculated as a linear combination of these induced 
currents. Next, the singular value decomposition is carried out 
and only those basis functions with a normalized singular 
value above T 1074 are retained. In other words, this 
procedure is equivalent to calculating the MBFs for a problem 
with a single supporting domain so that any potential source 
of error due to domain extensions is avoided. 

In order to choose the number of plane waves, the number 
of incident angles along 6 is set to N = |ka| + 10, where a 
is the radius of the minimum sphere enclosing the geometry. 
It is important to note that this choice, which is similar to the 
conventional rule for spherical wave truncation [70], fulfills 
the Nyquist criterion [69]. A similar discretization of the plane 
wave spectrum is accomplished for each azimuthal circle for 
a given 0. 

The number of surviving MBFs for different spheres and a 
cylinder is shown in Fig. 10. For the cylinder case, the height 
is set equal to the radius and no caps are considered (see inset 
in Fig. 10). In order to check the accuracy of the generated 
MBEFs, the current induced by a linearly polarized plane wave 
is compared with the current computed using the MoM. For 
the cylinder case, the plane-wave is assumed to be polarized 
parallel to the axis of the cylinder and the propagation vector 
is orthogonal to the aforementioned axis. The error is defined 
as: 


|x aE _ xon | 
2 





(20) 


e= - 
poeg 


MoM MBF 


where x and x are the coefficients of the low-level 
basis functions using the MoM and the single-block MBF 
approach. For the sphere case, the error ranges from 5.6- 1074 
to 7.7- 10-4 whereas this error ranges from 1.3 - 1074 to 
10-1074 for the cylinder. This verifies that the generated MBFs 
are good for modeling an arbitrary induced current. 
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Fig. 10. DoFs for several canonical geometries. 


The results presented in Fig. 10 reveal that the asymptotic 
behavior for large domains is approximately linear. It is 
interesting to note that the discretization of a geometry with 
low-level basis functions also involved a linear increase in 
the number of unknowns versus the surface to be analyzed. 
However, the slopes of the plots with MBFs are considerably 
smaller than when we consider low-level basis functions 
whose proportion is on the order of 100 basis functions per 
square wavelength. For instance, a sphere with a radius of 
2X (i.e., surface 50.27\7) can be modeled with 1056 MBFs 
whereas the number of low-level basis functions required 
would be on the order of 5027. 

Another interesting fact inferred from Fig. 10 pertains to 
the number of MBFs needed for electrically small domains. 
In this case, the asymptotic behavior is not reached and the 
number of MBFs grows much faster than in the case of 
electrically larger surfaces. In order to illustrate the validity 
of the latter observation, let us consider the sphere geometry 
once again. Furthermore, the relationship between the number 
of MBFs and domain surface for any arbitrary surface is 
considered to be equal to that for the sphere. This assumption 
is reasonable since the behavior of the curves in Fig. 10 
has also been observed for other geometries such as plates 
or cubes. Under the previous hypotheses, an MBF method 
working with domains of 28.9)? will require approximately 
2.47 times fewer unknowns than needed using the same MBFs 
method working with domains of 3.27. However, increasing 
the size of the domain surface requires the handling of much 
larger blocks and, consequently, the computationally burden 
is significantly increased in this case. In the next section, 
techniques to mitigate this computational burden are described. 

Finally, it is worthwhile to remark that similar observations 
have also been made by a number of different authors. For ex- 
ample, this fact was observed in [71] and exploited to generate 
MBFs on large blocks by applying physical optics. Similarly, 
the analysis of some particular geometries for different block 
sizes revealed that fewer MBFs are required when employing 
larger blocks to achieve the same prescribed error level [33]. 
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B. Multilevel MBF approach 


The aforementioned behavior of the DoFs versus the domain 
size suggests that the best compression rates are achieved for 
electrically large subdomains. Nevertheless, generating MBFs 
for these subdomains can be computationally expensive. More- 
over, the generation becomes more demanding when applying 
the plane wave spectrum approach since the currents on each 
subdomain has to be solved for multiple right hand sides. In 
order to benefit from the reduction in number of unknowns, 
while maintaining the computational burden at a reasonable 
level, a multilevel scheme which is referred to as Multilevel 
Characteristic Basis Function Method (MLCBFM) [21], [22] 
has been proposed. This technique is based on a recursive gen- 
eration of the macrobasis functions. A hierarchical partitioning 
of the geometry is required before applying the multilevel 
approach. This step is illustrated in Fig. 11 for the NASA 
almond in the context of a two-level scheme. In this figure, 
a gap has been introduced between geometry partitions to 
emphasize the domains. Next, the macrobasis functions are 
generated from the bottom to the top level. At each level, 
the basis functions are expressed as linear combinations of 
the basis functions defined on the level underneath. Let us 
consider a multilevel MBF at the /-th level, which is denoted 
by F. The number of subdomains of level / — 1 inside the 
domain is given by N while the number of MBFs for the n 
subdomain is given by C (n). Then, according to the multilevel 
definition of the MBFs, F can be expressed as: 


N C(n) 
oD Fa 21) 
n=1 i=1 
where is are the weights of the linear combinations that 


are computed in the process of generation of the MBFs. 

Once the coefficients of the MBFs have been computed at 
each level, the matrix of the system of equations and the right 
hand side can be recursively computed from bottom to top by 
carrying out pre- and post-multiplications by the coefficients 
of the matrices containing the coefficients of the MBFs at the 
corresponding level. 

The above formulation can be applied to an arbitrary num- 
ber of levels to benefit from the compression rate associated 
with large blocks. Nevertheless, the number of MBFs for 
very large blocks, which is assumed to be independent of 
the number of underlying levels, prevents us from handling 
extremely large blocks. In other words, as long as the MBFs 
are correctly generated at each level, the number of DoFs 
at the top level is expected to depend only on the domain 
extension. Consequently, very large domain extensions involve 
a large number of DoFs and, therefore, upper levels cannot be 
efficiently handled because the number of underlying MBFs 
becomes very large. 

Our experience is that a two-level scheme provides a good 
trade-off between accuracy and compression rate without sig- 
nificantly increasing the computational burden. There appears 
to be a consensus among other authors [21], [22], [72], [73] 
on this issue, since they also limited the implementation to 
two levels. 
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Fig. 11. Multilevel partitioning of the NASA almond. 


In order to illustrate the multilevel approach, the bistatic 
radar cross section (RCS) of a helicopter-type geometry is 
considered. The problem is discretized at 400 MHz by using 
96844 RWG basis functions. Next, the surface is partitioned so 
that the low-level basis functions are grouped into 420 domains 
for the first-level (i.e., conventional) MBFs (see Fig. 12). 
Generating the corresponding MBFs by using the plane-wave 
spectrum approach, followed by an SVD thresholding with 
7 = 1074, results in a problem with 13275 unknowns. Next, 
these first-level MBFs are grouped into three second-level 
domains which corresponds to splitting the geometry into three 
blocks along the x-axis. As a consequence, the final number of 
unknowns (i.e., second-level MBFs) is reduced to only 3192. 

The bistatic radar cross section for the case of an incident 
field E; = exp(j kx) is shown in Fig. 12. The results 
obtained using the MLFMA implemented in the commercial 
software FEKO [74] are also shown for comparison purposes. 
The agreement between the two methods is seen to be excel- 
lent. 


VIII. CONCLUSIONS AND OUTLOOK 


Integral-equation approaches remain among the most com- 
petitive methods for the solution of large radiation or scattering 
problems. Domain decomposition started to be applied to 
this type of methods about fifteen years ago and has been 
introduced at about the same time by different labs, in several 
variants. They rely on the a priori determination of the 
subspace in which current distributions (or equivalent currents) 
on a given subdomain can be found. Those subspaces are 
subtended by Macro Basis Functions (MBFs), defined in terms 
of the original elementary basis functions. Their main advan- 
tage is that they allow the direct solution of large problems, 
thereby avoiding the uncertainty about the number of iterations 
and enabling very quick solutions for multiple excitations. 
We reviewed different techniques for the determination of the 
MBFs, as well as methods for the very fast calculation of 
their reactions. We also reviewed the similarities between this 
class of methods and Krylov-based iterative techniques, like 
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Fig. 12. Bistatic RCS for the helicopter. The inset shows the first level 
partitioning. 


the Full Orthogonalization Method (FOM). Besides the direct- 
solution aspect, the domain-decomposition nature of MBF- 
based methods allows one to enforce boundary-type conditions 
on each subdomain individually. 

This type of methods is particularly well-suited for appli- 
cation to geometries for which the characteristic dimensions 
of the mesh are very much finer than the wavelength, such 
that the number of physical degrees of freedom of the fields 
is much smaller than those implied by the complexity of 
the geometrical discretization. Such geometries are extremely 
common in industrial applications. MBFs can be applied to 
preconditioned systems of equations, which can implicitly 
account for the connectivity between subdomains without 
significant loss of accuracy. Besides, beyond this precondi- 
tioning, the MBF approach in itself seems to add its own 
preconditioning effect, such that practically no loss of accuracy 
is observed when solving scattering problems on objects near 
resonance. Being non-iterative in nature, the MBF approach 
lends itself very well to parallelization, which means that 
the MBF methodology is among the approaches that will at 
best benefit from the current strong trend toward large-scale 
multiple-core calculation. 

It is expected that improved construction of MBFs, faster 
interaction methods and multi-level approaches will receive 
more attention in the coming years. As previously pointed 
out, the generation of MBFs yields a different number of 
unknowns, depending on the employed scheme. Problems 
requiring a restricted set of excitations can benefit from the 
high compression provided by the primary-secondary approach 
(or its extension to higher orders) without a significant loss of 
accuracy. On the other side, the plane-wave approach enables 
us to compute excitation-free MBFs at the cost of an increment 
in the number of degrees of freedom. This fact suggests that 
a priori knowledge of the excitation can be employed to 
decrease the degrees of freedom in the problem. 

Moreover, considering larger domains in the MBF genera- 
tion seems to (proportionally) require fewer unknowns, setting 
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the basis for the multilevel formulation. This probably results 
from the fact that more information is used during the MBFs 
generation and, consequently, the achieved MBFs can better 
model the final currents. Thus, a rigorous study is still needed 
regarding the minimum number of degrees of freedom to be 
associated with a given problem, and regarding how to fully 
exploit all the a priori information about optimum domain size 
and excitation. Regarding the latter, for instance, excitation- 
dependent MBFs may remain valid when the primary source 
is moved over a finite domain. MBFs may also have a regular- 
izing effect in the solution of radiation or scatering problems. 
Besides providing more stable solutions, their physics-based 
foundation can also serve purposes that are other than purely 
computational. For instance, they have been exploited for 
calibration purposes in [75]-[78]. They could also be exploited 
with advantage in time-domain solution schemes, since MBF- 
type basis functions in the time domain have been observed 
to provide a better late-time stability [79]. 

To conclude, MBF methods and related techniques have 
already proved quite powerful for the solution of Maxwell’s 
equations in surface integral-equation form. Over the past 
decade, they have been strongly accelerated and methodologies 
to produce more complete sets of MBFs have been developed. 
Their direct-solution nature and their ease of parallelization 
make them preferable to iterative techniques for a wide class of 
problems. Nevertheless, it is expected that combined research 
on MBF-based and iterative techniques will be beneficial to the 
development of both solution methodologies. Further progress 
in this area may also benefit from continued research on the 
physical degrees of freedom of fields excited on arbitrarily- 
shaped structures, with possible applications beyond those in 
computational methods. 
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